{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Copyright 2020 The OpenFermion Developers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#@title Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "# you may not use this file except in compliance with the License.\n",
    "# You may obtain a copy of the License at\n",
    "#\n",
    "# https://www.apache.org/licenses/LICENSE-2.0\n",
    "#\n",
    "# Unless required by applicable law or agreed to in writing, software\n",
    "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "# See the License for the specific language governing permissions and\n",
    "# limitations under the License."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hamiltonian Time Evolution and Expectation Value Computation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table class=\"tfo-notebook-buttons\" align=\"left\">\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://quantumai.google/openfermion/fqe/tutorials/hamiltonian_time_evolution_and_expectation_estimation\"><img src=\"https://quantumai.google/site-assets/images/buttons/quantumai_logo_1x.png\" />View on QuantumAI</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://colab.research.google.com/github/quantumlib/OpenFermion/blob/master/docs/fqe/tutorials/hamiltonian_time_evolution_and_expectation_estimation.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/colab_logo_1x.png\" />Run in Google Colab</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://github.com/quantumlib/OpenFermion/blob/master/docs/fqe/tutorials/hamiltonian_time_evolution_and_expectation_estimation.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/github_logo_1x.png\" />View source on GitHub</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a href=\"https://storage.googleapis.com/tensorflow_docs/OpenFermion/docs/fqe/tutorials/hamiltonian_time_evolution_and_expectation_estimation.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/download_icon_1x.png\" />Download notebook</a>\n",
    "  </td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial describes the FQE's capabilities for Hamiltonian time-evolution and expectation value estimation\n",
    "\n",
    "Where possible, LiH will be used as an example molecule for the API."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    import fqe\n",
    "except ImportError:\n",
    "    !pip install fqe --quiet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Print = True\n",
    "from openfermion import FermionOperator, MolecularData\n",
    "from openfermion.utils import hermitian_conjugated\n",
    "import numpy\n",
    "import fqe\n",
    "from fqe.unittest_data import build_lih_data, build_hamiltonian\n",
    "numpy.set_printoptions(floatmode='fixed', precision=6, linewidth=80, suppress=True)\n",
    "numpy.random.seed(seed=409)\n",
    "\n",
    "h1e, h2e, wfn = build_lih_data.build_lih_data('energy')\n",
    "lih_hamiltonian = fqe.get_restricted_hamiltonian(([h1e, h2e]))\n",
    "lihwfn = fqe.Wavefunction([[4, 0, 6]])\n",
    "lihwfn.set_wfn(strategy='from_data', raw_data={(4, 0): wfn})\n",
    "if Print:\n",
    "    lihwfn.print_wfn()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application of one- and two-body fermionic gates\n",
    "\n",
    "The API for time propogation can be invoked through the fqe namespace or the wavefunction object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dummy geometry\n",
    "from openfermion.chem.molecular_data import spinorb_from_spatial\n",
    "from openfermion import jordan_wigner, get_sparse_operator, InteractionOperator, get_fermion_operator\n",
    "\n",
    "h1s, h2s = spinorb_from_spatial(h1e, numpy.einsum(\"ijlk\", -2 * h2e) * 0.5)\n",
    "mol = InteractionOperator(0, h1s, h2s)\n",
    "ham_fop = get_fermion_operator(mol)\n",
    "ham_mat = get_sparse_operator(jordan_wigner(ham_fop)).toarray()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from scipy.linalg import expm\n",
    "time = 0.01\n",
    "evolved1 = lihwfn.time_evolve(time, lih_hamiltonian)\n",
    "if Print:\n",
    "    evolved1.print_wfn()\n",
    "evolved2 = fqe.time_evolve(lihwfn, time, lih_hamiltonian)\n",
    "if Print:\n",
    "    evolved2.print_wfn()\n",
    "assert numpy.isclose(fqe.vdot(evolved1, evolved2), 1)\n",
    "cirq_wf = fqe.to_cirq_ncr(lihwfn)\n",
    "evolve_cirq = expm(-1j * time * ham_mat) @ cirq_wf\n",
    "test_evolve = fqe.from_cirq(evolve_cirq, thresh=1.0E-12)\n",
    "assert numpy.isclose(fqe.vdot(test_evolve, evolved1), 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exact evolution implementation of quadratic Hamiltonians\n",
    "\n",
    "Listed here are examples of evolving the special Hamiltonians.\n",
    "\n",
    "Diagonal Hamiltonian evolution is supported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "wfn = fqe.Wavefunction([[4, 2, 4]])\n",
    "wfn.set_wfn(strategy='random')\n",
    "if Print:\n",
    "    wfn.print_wfn()\n",
    "\n",
    "diagonal = FermionOperator('0^ 0', -2.0) + \\\n",
    "           FermionOperator('1^ 1', -1.7) + \\\n",
    "           FermionOperator('2^ 2', -0.7) + \\\n",
    "           FermionOperator('3^ 3', -0.55) + \\\n",
    "           FermionOperator('4^ 4', -0.1) + \\\n",
    "           FermionOperator('5^ 5', -0.06) + \\\n",
    "           FermionOperator('6^ 6', 0.5) + \\\n",
    "           FermionOperator('7^ 7', 0.3)\n",
    "if Print:\n",
    "    print(diagonal)\n",
    "    \n",
    "evolved = wfn.time_evolve(time, diagonal)\n",
    "if Print:\n",
    "    evolved.print_wfn()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exact evolution of dense quadratic hamiltonians is supported.  Here is an evolution example using a spin restricted Hamiltonian on a number and spin conserving wavefunction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "norb = 4 \n",
    "h1e = numpy.zeros((norb, norb), dtype=numpy.complex128) \n",
    "for i in range(norb): \n",
    "    for j in range(norb): \n",
    "        h1e[i, j] += (i+j) * 0.02 \n",
    "    h1e[i, i] += i * 2.0 \n",
    "\n",
    "hamil = fqe.get_restricted_hamiltonian((h1e,)) \n",
    "wfn = fqe.Wavefunction([[4, 0, norb]]) \n",
    "wfn.set_wfn(strategy='random') \n",
    "initial_energy = wfn.expectationValue(hamil) \n",
    "print('Initial Energy: {}'.format(initial_energy))\n",
    "evolved = wfn.time_evolve(time, hamil) \n",
    "final_energy = evolved.expectationValue(hamil)\n",
    "print('Final Energy:   {}'.format(final_energy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The GSO Hamiltonian is for evolution of quadratic hamiltonians that are spin broken and number conserving."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "norb = 4 \n",
    "h1e = numpy.zeros((2*norb, 2*norb), dtype=numpy.complex128) \n",
    "for i in range(2*norb): \n",
    "    for j in range(2*norb): \n",
    "        h1e[i, j] += (i+j) * 0.02 \n",
    "    h1e[i, i] += i * 2.0 \n",
    "\n",
    "hamil = fqe.get_gso_hamiltonian((h1e,)) \n",
    "wfn = fqe.get_number_conserving_wavefunction(4, norb) \n",
    "wfn.set_wfn(strategy='random') \n",
    "initial_energy = wfn.expectationValue(hamil) \n",
    "print('Initial Energy: {}'.format(initial_energy))\n",
    "evolved = wfn.time_evolve(time, hamil) \n",
    "final_energy = evolved.expectationValue(hamil)\n",
    "print('Final Energy:   {}'.format(final_energy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The BCS hamiltonian evovles spin conserved and number broken wavefunctions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "norb = 4\n",
    "time = 0.001\n",
    "wfn_spin = fqe.get_spin_conserving_wavefunction(2, norb)\n",
    "hamil = FermionOperator('', 6.0)\n",
    "for i in range(0, 2*norb, 2):\n",
    "    for j in range(0, 2*norb, 2):\n",
    "        opstring = str(i) + ' ' + str(j + 1)\n",
    "        hamil += FermionOperator(opstring, (i+1 + j*2)*0.1 - (i+1 + 2*(j + 1))*0.1j)\n",
    "        opstring = str(i) + '^ ' + str(j + 1) + '^ '\n",
    "        hamil += FermionOperator(opstring, (i+1 + j)*0.1 + (i+1 + j)*0.1j)\n",
    "h_noncon = (hamil + hermitian_conjugated(hamil))/2.0\n",
    "if Print:\n",
    "    print(h_noncon)\n",
    "\n",
    "wfn_spin.set_wfn(strategy='random')\n",
    "if Print:\n",
    "    wfn_spin.print_wfn()\n",
    "\n",
    "spin_evolved = wfn_spin.time_evolve(time, h_noncon)\n",
    "if Print:\n",
    "    spin_evolved.print_wfn()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exact Evolution Implementation of Diagonal Coulomb terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "norb = 4\n",
    "wfn = fqe.Wavefunction([[5, 1, norb]])\n",
    "vij = numpy.zeros((norb, norb, norb, norb), dtype=numpy.complex128)\n",
    "for i in range(norb):\n",
    "            for j in range(norb):\n",
    "                vij[i, j] += 4*(i % norb + 1)*(j % norb + 1)*0.21\n",
    "                \n",
    "wfn.set_wfn(strategy='random')\n",
    "\n",
    "if Print:\n",
    "    wfn.print_wfn()\n",
    "    \n",
    "hamil = fqe.get_diagonalcoulomb_hamiltonian(vij)\n",
    "    \n",
    "evolved = wfn.time_evolve(time, hamil)\n",
    "if Print:\n",
    "    evolved.print_wfn()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exact evolution of individual n-body anti-Hermitian gnerators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "norb = 3\n",
    "nele = 4\n",
    "ops = FermionOperator('5^ 1^ 2 0', 3.0 - 1.j)\n",
    "ops += FermionOperator('0^ 2^ 1 5', 3.0 + 1.j)\n",
    "wfn = fqe.get_number_conserving_wavefunction(nele, norb)\n",
    "wfn.set_wfn(strategy='random')\n",
    "wfn.normalize()\n",
    "if Print:\n",
    "    wfn.print_wfn()\n",
    "evolved = wfn.time_evolve(time, ops)\n",
    "if Print:\n",
    "    evolved.print_wfn()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Approximate evolution of sums of n-body generators\n",
    "\n",
    "Approximate evolution can be done for dense operators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "lih_evolved = lihwfn.apply_generated_unitary(time, 'taylor', lih_hamiltonian, accuracy=1.e-8)\n",
    "if Print:\n",
    "    lih_evolved.print_wfn()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "norb = 2\n",
    "nalpha = 1\n",
    "nbeta = 1\n",
    "nele = nalpha + nbeta\n",
    "time = 0.05\n",
    "h1e = numpy.zeros((norb*2, norb*2), dtype=numpy.complex128)\n",
    "for i in range(2*norb):\n",
    "    for j in range(2*norb):\n",
    "        h1e[i, j] += (i+j) * 0.02\n",
    "    h1e[i, i] += i * 2.0\n",
    "hamil = fqe.get_general_hamiltonian((h1e,))\n",
    "spec_lim = [-1.13199078e-03, 6.12720338e+00]\n",
    "wfn = fqe.Wavefunction([[nele, nalpha - nbeta, norb]])\n",
    "wfn.set_wfn(strategy='random')\n",
    "if Print:\n",
    "    wfn.print_wfn()\n",
    "evol_wfn = wfn.apply_generated_unitary(time, 'chebyshev', hamil, spec_lim=spec_lim)\n",
    "if Print:\n",
    "    evol_wfn.print_wfn()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "API for determining desired expectation values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "rdm1 = lihwfn.expectationValue('i^ j')\n",
    "if Print:\n",
    "    print(rdm1)\n",
    "val = lihwfn.expectationValue('5^ 3')\n",
    "if Print:\n",
    "    print(2.*val)\n",
    "trdm1 = fqe.expectationValue(lih_evolved, 'i j^', lihwfn)\n",
    "if Print:\n",
    "    print(trdm1)\n",
    "val = fqe.expectationValue(lih_evolved, '5 3^', lihwfn)\n",
    "if Print:\n",
    "    print(2*val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2.B.1 RDMs \n",
    "In addition to the above API higher order density matrices in addition to hole densities can be calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "rdm2 = lihwfn.expectationValue('i^ j k l^')\n",
    "if Print:\n",
    "    print(rdm2)\n",
    "rdm2 = fqe.expectationValue(lihwfn, 'i^ j^ k l', lihwfn)\n",
    "if Print:\n",
    "    print(rdm2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2.B.2 Hamiltonian expectations (or any expectation values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "li_h_energy = lihwfn.expectationValue(lih_hamiltonian)\n",
    "if Print:\n",
    "    print(li_h_energy)\n",
    "li_h_energy = fqe.expectationValue(lihwfn, lih_hamiltonian, lihwfn)\n",
    "if Print:\n",
    "    print(li_h_energy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2.B.3 Symmetry operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "op = fqe.get_s2_operator()\n",
    "print(lihwfn.expectationValue(op))\n",
    "op = fqe.get_sz_operator()\n",
    "print(lihwfn.expectationValue(op))\n",
    "op = fqe.get_time_reversal_operator()\n",
    "print(lihwfn.expectationValue(op))\n",
    "op = fqe.get_number_operator()\n",
    "print(lihwfn.expectationValue(op))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
